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Abstract.  This  paper  is  concerned  with  a  stochastic  inverse  methodology  arising  in 
electromagnetic  imaging.  Nondestructive  testing  using  guided  microwaves  covers  a  wide 
range  of  industrial  applications  including  early  detection  of  anomalies  in  conducting  ma¬ 
terials.  The  focus  of  this  paper  is  the  identification  of  electromagnetic  material  parame¬ 
ters  and  emphasis  is  on  one  dimensional  scattering  of  a  dielectric  slab.  The  direct  problem 
can  be  solved  numerically  using  the  finite-difference  time-domain  method  (FDTD).  The 
Markov  Chain  Monte  Carlo  method  (MCMC)  is  applied  to  the  inversion  problem.  Some 
successful  results  of  computational  experiments  are  demonstrated  in  order  to  show  the 
feasibility  and  applicability  of  the  proposed  method. 

Keywords:  Microwave,  Dielectric  loss,  Electrical  cable,  FDTD,  Metropolis-Hasting  al¬ 
gorithm 


1.  Introduction.  Recently,  interest  has  been  growing  in  structural  health  monitoring 
(SHM)  related  to  aging  management  of  large  scale  systems,  such  as  airplanes,  bridges 
and  nuclear  power  plants.  Various  kinds  of  nondestructive  testing  (NDT)  techniques  such 
as  ultrasonic,  eddy  current,  thermal  are  applied  to  the  detection  and  the  characteriza¬ 
tion  of  material  damage.  Combining  NDT  with  simulation  is  a  key  component  of  future 
structural  monitoring  technologies.  Mathematical  descriptions  of  non-destructive  test¬ 
ing  (NDT)  can  be  formulated  as  either  a  forward  or  an  inverse  problem  for  the  physical 
domain  of  the  inspection.  A  forward  problem  represents  a  real  NDT  system  mathemat¬ 
ically  using  the  input  and  the  output  relationship  with  the  appropriate  admissible  class 
of  material  parameters,  while  an  inverse  problem  involves  the  construction  of  a  method 
for  recovering  and/or  visualizing  material  flaw  information  using  the  mathematical  for¬ 
mulation  of  the  forward  problem.  Figure  1  demonstrates  the  overall  configuration  of  the 
proposed  system.  In  this  paper,  a  stochastic  inverse  methodology  for  NDT  arising  in 
electromagnetic  imaging  is  investigated.  Nondestructive  testing  using  guided  microwaves 
is  used  in  a  wide  range  of  industrial  applications  including  early  detection  of  anomalies  in 
conducting  materials,  analysis  of  human  muscle  tissue,  etc.  The  characterization  of  mate¬ 
rial  corrosion  damage  continues  to  be  a  very  challenging  problem.  The  problem  is  further 
complicated  by  the  dispersive  nature  of  the  insulation  covering  the  corrosion.  Estimating 
the  material  properties  in  conventional  measurements  has  been  studied  extensively  [1]. 
For  instance,  L.  F.  Chen  et  al.  have  addressed  the  practical  guidance  on  the  development 
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Scanning  Strategies 


Figure  1.  Diagnosis  system  using  forward  and  inverse  analysis 


of  suitable  measurement  methodologies  using  micro  guided  wave  for  a  variety  of  materi¬ 
als  in  their  book  [1],  Our  focus  in  this  paper  is  in  the  identification  of  electromagnetic 
material  parameters  and  the  emphasis  is  on  one  dimensional  scattering  of  a  dielectric 
slab.  Although  prior  work  exists  using  nonlinear  least  square  methods  [2,  3,  4,  5,  6],  it  is 
well  known  that  the  problem  mentioned  above  has  many  solutions  due  to  the  fact  that 
it  is  ill-posed.  Recently,  interest  has  grown  in  stochastic  inversion  using  Markov  Chain 
Monte  Carlo  (MCMC)  methods  [7,  8].  Although  the  book  written  by  L.  F.  Chen  [1]  et 
al.  contains  around  thousands  of  references  in  a  variety  of  applications,  there  is  not  any 
literature  using  MCMC  to  characterize  the  material  properties.  The  author  first  proposed 
the  parameter  estimation  method  using  the  Gibbs  sampling  algorithm  in  [9] .  The  method 
has  great  advantages  for  the  more  practical  aspects  of  inversions  such  as  setting  initial 
guesses  and  overcoming  local  minimums  of  the  inverse  solution.  In  this  paper,  a  stochastic 
inversion  technique  based  on  the  Metropolis-Hasting  algorithm  is  applied  to  our  problem. 


2.  Formulation  of  the  Inverse  Problem.  Let  Dx(t,z),  Hy(t,z )  and  Ex(t,z )  be  the 
component  of  the  electric  flux  density,  and  the  magnetic  and  electric  field  at  time  t  G  [0,  T] 
and  at  location  z  e  [0,  Z\.  Using  normalized  electric  flux  densities  and  electric  fields 


D,r  = 


\AoFo 


-.D, 


(1) 


and  from  Maxwell  equations,  the  electromagnetic  propagation  in  ^-direction  is  governed 

by 


dDx(t,z ) 
dt 

Dx(uj,z) 
dHy(t,z ) 
dt 


1  3Hy(t,Z ) 

uVeoAD  dz 

e*r(uj)Ex(uj,z ) 

1  dEx(t,z ) 

yToFo  dz 


+  Js(t ,  z) 


in  [0,  T]  x 


in  [0,  T\  x  [0,  Z] 


M 


(2) 

(3) 

(4) 
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Figure  2.  Dielectric  constant  for  a  Debye  medium  over  the  frequency 
range  of  10  to  1000 [MHz\ 


with  zero  initial  states  and  absorbing  boundary  conditions.  In  system  (2)  and  (4),  eo  and 
Ho  denote  the  dielectric  constant  and  magnetic  permeability  of  vacuum,  respectively. 

The  source  current  Js  is  the  test  signal  in  the  inspection  process  and  is  given  by 


Js(t,  z )  =  5(z  -  z8)gs{t)I(p,t.){t) 


(5) 


where  5  denote  the  Dirac  distribution  and  zs  G  [0,  Z\)  denotes  the  source  point  of  test 
signal.  The  test  signal  is  truncated  at  a  finite  time  ts  by  the  indicator  function  /. 

The  spatial  domain  [0,  Z\  is  decomposed  into  two  separate  regions.  The  region  z  G 
[0 ,Zi)  is  air  and  is  assumed  to  have  zero  electric  polarization  and  zero  conductivity. 
Hence,  in  (3),  it  becomes 

<  =  I-  (6) 


A  target  material  in  z  G  [Z\,  Z]  is  assumed  to  be  a  chemical  material  such  as  a  polymer 
where  the  dielectric  constant  and  the  conductivity  are  dispersive.  A  material  like  this  is 
well  approximated  by  the  Debye  law  described  by 


e*  un)  —  er  + 


a 


+ 


Xi 


juje0  1  +  jujT 


(7) 


where  er,  a,  xi  and  t  denote  a  dielectric  constant,  conductivity,  and  frequency-dependent 
parameters  [3].  As  indicated  in  Equation  (7),  a  target  material  includes  a  medium  whose 
dielectric  constant  and  conductivity  vary  over  the  frequency  range.  Figures  2  and  3 
show  a  relative  dielectric  constant  and  conductivity  as  functions  of  frequency  for  a  Debye 
medium  (7)  with  properties  (21).  From  (7),  it  follows  that  the  dielectric  constant  for  a 
Debye  medium  can  be  represented  as 


er  + 


Xi 

1  +  (tcj)2 


and  the  frequency  dependent  conductivity  can  be  rewritten  as 


respectively. 


(ne0 


cr 

e0u>  1  +  (tcu)2 


> 
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Figure  3.  Conductivity  for  a  Debye  medium  over  the  frequency  range  of 
10  to  1000  [MHz] 


There  are  many  interesting  applications  for  identifying  frequency  dependent  parameters 
in  the  Debye  model  (7).  Thus  our  inverse  problem  is  to  identify  dielectric  parameters 

q  =  {er,cr,xi,T} 

from  observations  near  the  frequency  dependent  medium: 

Y (£;  q)  =  Ex(t,  0;  Js)  (8) 

where  Ex(t,z)  is  the  solution  of  (2)- (4). 


3.  Numerical  Scheme  of  Direct  Problem.  The  detection  signal  (8)  corresponding 
to  (5)  can  be  obtained  by  numerical  simulation  using  the  finite-difference  time-domain 
(FDTD)  method  [10].  FDTD  uses  central-difference  approximations  to  the  space  and 
time  partial  derivatives.  The  resulting  finite-difference  equations  are  solved  in  a  so-called 
“leapfrog  manner” .  More  precisely,  the  electrical  field  components  in  a  volume  of  space 
are  solved  at  a  given  instant  in  time.  Then  the  magnetic  hied  components  in  the  same 
spatial  volume  are  solved  at  the  next  instant  in  time. 

Taking  the  central  difference  formula  for  both  the  temporal  and  spatial  derivatives, 
Equations  (2)  and  (4)  can  be  approximated  by 


Pnx+l,2[k\- bnx-1/2[k] 

At 

_ 1  H;lk  +  1/2]  -  HJlk  -  1/2]  + 

H£+l[k  +  1/2]  -  H£[k  +  1/2] 

At 

1  Ex+1/2[k  +  l}~  Ex+1/2[k\ 

Veoho  Ar 


(9) 

(10) 


where  n  implies  time  division,  i.e. ,  t  =  At- n  and  where  the  terms  in  parentheses  represent 
grid  point  in  the  spatial  domain,  i.e.,  z  =  Az-k,  respectively.  In  a  one  dimensional  spatial 
domain,  there  exist  two  boundary  points  z  —  0,  Z.  An  absorption  boundary  condition 
is  necessary  at  both  end  points,  z  =  0,  Z  in  order  to  keep  the  outgoing  electrical  and 
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magnetic  fields  E  and  H  from  being  reflected  back  into  the  target  domain.  This  can  be 
accomplished  by  setting 

b;[0]  =  £>r‘[i]  (11) 

4"[o]  =  4”_1[i]  (12) 

bl\K]  =  b:~'{K-l]  (13) 

El\K\  =  ET'lK-l]  (14) 

where  Z  =  K  ■  A z.  To  transform  Equation  (3)  into  a  time  domain  difference  equation  for 

implementation  into  FDTD,  the  frequency  term  should  be  replaced  by  the  time  domain 
representation.  Noting  that  1  / juj  in  the  frequency  domain  is  equivalent  to  integration  in 
the  time  domain,  Equation  (3)  becomes 


a 


Dx  (t)  —  erEx(t)  H - / 

eo  Jo 


Er  ( t 


dt+ f 


exp 


t  -t 


ET  ( t  dt  . 


(15) 


By  approximating  Equation  (15)  as  a  summation  in  the  sampled  time  domain,  the  nu¬ 
merical  scheme  can  be  represented  by 

Dx  =  er  •  En  +  In  +  Sn  (16) 


where 


r  =  i 


n—  1 


+ 


cr  •  At 


■  El 


eo 

Sn  =  exp  ( 


sa-'  +  xiA-K 


(17) 

(18) 


starting  with  zero  initial  states  with  respect  to  both  D®  and  E®.  Consequently,  the  forward 
problem  can  be  implemented  by  the  repeated  scheme  mentioned  above: 


y"(  q)  =  J-K{0](j,)- 

V  eo 


(19) 


Figure  4  depicts  a  simulation  result  for  a  direct  problem  (19).  A  sinusoidal  wave  was 
applied  to  the  test  signal  at  zs  =  5  •  in  Equation  (5): 

9s  =  sin(2  •  7T  ■  /  •  At  ■  n)  (20) 

where  /  determines  the  frequency  of  the  test  signal.  In  the  experiment,  a  sinusoidal  wave 
of  700 MHz  is  applied  to  a  dielectric  medium  with  parameters 

er  =  2.0,  cr  =  0.01,  Xi  —  2-0,  r  =  0.001[/isec].  (21) 


4.  Stochastic  Inverse  Methodology.  The  output  least  square  method  is  a  conven¬ 
tional  inverse  methodology  which  seeks  the  optimal  solution  of 

K 

J (q)  =  minqeQ  lF(^  <l)  ~  Yd  f  •  (22) 

i=  1 

However,  it  is  well  known  that  those  experimental  tests  might  not  achieve  good  conver¬ 
gence  results  because  of  the  complicated  dependence  in  the  related  direct  problem  (8).  In 
this  paper,  we  propose  a  new  stochastic  inverse  methodology  for  identifying  those  parame¬ 
ters.  Using  a  Bayesian  formula,  the  full  probability  model  is  specified  by  the  deterministic 
Formula  (8).  The  posteriori  density  function  with  respect  to  the  set  of  unknown  parameter 
vector  can  be  represented  by 

vr(q|Yd)  oc  /(Yd|q)7r(q).  (23) 
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FIGURE  4.  Simulation  result 


Suppose  that  measurements  are  made 

Yd  =  Y  (U,  q)  +  Vi,  0,£2),  [i  =  1,2, . 

Then  the  likelihood  functional  can  be  written  by  the  following  form: 


K 


KYd\q)  oc  n  — ^=-  ex^ 

i=l  *  ’ 


|rfeq)-rjl 

2p 


(24) 


(25) 


In  our  estimation  algorithm,  prior  distributions  for  unknown  parameters  vr(q)  are  specified 
as 

4 

^(q)=n^)-  (2e) 

1=1 

Then  we  specify  full  conditional  distributions  for  posteriori  density  function  through  the 
likelihood  functionals.  The  full  conditional  distribution  Ti(qi\)  ( l  =  1,2, 3, 4)  has  the 
representation 


,  ,  ,n  1  (  \Y(t,;q)-Y‘ 

n{qi\q-hYd)  oc  7r(g) )_[_[  — exp  I - — - 

where  denotes  the  remaining  component  except  <37. 

Our  estimation  algorithm  is  based  on  sampling  procedures  from  the  posteriori  distri¬ 
bution  from  which  a  sample  can  be  drawn  via  Markov  chains.  To  this  end,  a  transition 
kernel  p(q,  (ft)  must  be  constructed  in  such  a  way  that  the  posteriori  function  Equation 
(23)  is  the  equilibrium  distribution  of  the  chain.  A  simple  way  to  implement  this  is  to 
consider  the  reversible  chains  where  the  kernel  p  satisfies 


7r(q)p(q,  (ft)  =  7 r((ft)p((ft,  q).  (28) 

The  feature  of  the  Metropolis-Hasting  algorithm  is  that  this  reversible  kernel  p  is  con¬ 
structed  by 


P(V,A)=  /  k{q,(ft)a(q,(ft)d(ft  +  I(q  e  A)  {  1  -  /  k(q,  (ft)a{q,  (ft)  \  deft 


'A 


(29) 
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TABLE  1 .  Estimation  summary  using  MH  algorithm 


Quantities 

True  Values 

Estimated  Values 

€r 

2.00 

1.99 

login  o’ 

-2.00 

-1.95 

Xi 

2.00 

1.68 

logio  T 

-9.00 

-9.05 

for  any  subset  A  of  the  parameter  space.  In  the  above  expression,  the  acceptance  proba¬ 
bility 

/  / \  •  J\  n((j))k(faq)  1  /or.'t 

aM)  =  mm\l’n(q)k(q,4.)S  <30) 

plays  an  essential  role  in  MH  algorithms  [8].  There  exist  varieties  of  the  proposal  kernel 
functions  k  for  practical  implementations.  Thus  the  process  version  of  the  MH  algorithm 
is  given  by  the  following  steps: 

Step  1:  Initialize  the  iteration  counter  j  —  1  and  set  the  initial  guess  of  the  chain  q°. 
Step  2  :  Initialize  the  component  i  —  1. 

Step  3  :  Move  the  ?'th  component  of  the  parameter  vector  of  the  chain  to  a  new  value  fa 
generated  from  the  prescribed  transition  kernel  k-L  (qf  ^,fa ). 


Step  4  :  Calculate  the  acceptance  probability  of  the  move  ccj  \  1\  fa  j  given  by  Equa¬ 

tion  (30).  If  the  move  is  accepted,  update  the  chain  qf}  =  fa.  If  the  move  is  rejected,  set 
(j)  0—1) 

Qi  =  Qi  ■ 


Step  5:  Change  the  counter  from  i  to  i  +  1  and  return  to  Step  3  until  the  dimension  of 
parameter  vector,  i.e. ,  dim{ q)  =  4. 

Step  6:  Change  the  counter  from  j  to  j  +  1  and  return  to  Step  2  until  convergence  is 
reached. 


5.  Simulation  Experiments.  Throughout  our  numerical  experiments,  simulation  data 
are  generated  using  the  direct  problem  in  Section  3.  Taking  into  account  that  actual 
measurement  procedures  can  be  well  approximated  by  Equation  (24),  some  random  noise 
is  added  to  the  simulation  data  to  test  our  inverse  methodology.  The  second  variance  term 
in  Equation  (24)  was  provided  using  a  conventional  Gaussian  random  generator.  The  set 
of  true  parameters  in  the  experiments  is  selected  by  Equation  (21).  Table  1  shows  the 
estimated  results  in  the  experiments. 

6.  Conclusion.  The  inverse  methodology  arising  in  electromagnetic  imaging  was  dis¬ 
cussed  within  the  Bayesian  framework.  We  propose  the  statistical  method  for  identifying 
electromagnetic  material  parameters  on  a  one  dimensional  scattering  problem  of  a  dielec¬ 
tric  slab.  The  direct  problem  was  obtained  by  numerical  simulation  using  FDTD  method. 
The  stochastic  inverse  problem  was  solved  via  the  Metropolis-Hasting  algorithm. 

Acknowledgment.  This  material  is  based  on  research  sponsored  by  AOARD,  under 
agreement  number  No.  FA2386-10-1-4076. 
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